﻿using System;
using System.Collections.Generic;
using System.Drawing;
using System.Text;
using System.Windows.Forms;

namespace Mclib
{
	/// <summary>
	/// 
	/// </summary>
	public class IsoContours : IComparable
	{
		/// <summary>
		/// Radians of 1/4 circle.
		/// </summary>
		public readonly double RADIANS_1_4 = 0.5*Math.PI;
		/// <summary>
		/// Radians of 3/4 circle.
		/// </summary>
		public readonly double RADIANS_3_4 = 1.5*Math.PI;
		/// <summary>
		/// Radians of 1/8 circle.
		/// </summary>
		public readonly double RADIANS_1_8 = 0.25*Math.PI;
		/// <summary>
		/// Radians of 3/8 circle.
		/// </summary>
		public readonly double RADIANS_3_8 = 0.75*Math.PI;
		/// <summary>
		/// Radians of 5/8 circle.
		/// </summary>
		public readonly double RADIANS_5_8 = 1.25*Math.PI;
		/// <summary>
		/// Radians of 7/8 circle.
		/// </summary>
		public readonly double RADIANS_7_8 = 1.75*Math.PI;
		/// <summary>
		/// Square-root of one half.
		/// </summary>
		public readonly double ROOT_1_2 = Math.Sqrt(0.5);
		/// <summary>
		/// Square-root of two.
		/// </summary>
		public readonly double ROOT_2 = Math.Sqrt(2.0);
		/// <summary>
		/// The level of 'Z' for which these contours are defined.
		/// </summary>
		public double IsoZ;
		/// <summary>
		/// The data for which the contour set is defined.  Z[i,j] is a tabulation of the function Z(y,x).
		/// </summary>
		public double[,] Z;
		/// <summary>
		/// The complete list contours that trace an IsoZ value sorted from longest to shortest.
		/// </summary>
		public List<Contour> Contours = new List<Contour>();
		/// <summary>
		/// The vertex pairs computed as the contours were being built.
		/// </summary>
		public List<ContourSegment> Segments = new List<ContourSegment>();
		/// <summary>
		/// The vertices computed as the conoturs were being built.
		/// </summary>
		public List<ContourVertex> Vertices = new List<ContourVertex>();
		/// <summary>
		/// The rectangle that delimits the extent of the histogram mesh.
		/// </summary>
		public RectangleF MeshRect;
		/// <summary>
		/// Constructs the iso-value contours for the function tabulated in Z.
		/// </summary>
		/// <param name="Z"></param>
		/// <param name="IsoZ">The value of Z isolated by the contours.</param>
		public IsoContours(double[,] Z, double IsoZ)
		{
			this.Z = Z;
			this.IsoZ = IsoZ;
			int iLen = Z.GetLength(0);
			int jLen = Z.GetLength(1);
			MeshRect = RectangleF.FromLTRB(0.0f,0.0f,Z.GetLength(1)-1,Z.GetLength(0)-1);
			ContourVertex[,] zCrossings = new ContourVertex[iLen,jLen];
			double u0,v0,u1,v1,z00,z01,z10,z11,dz0,dz1;
			bool[] crosses = new bool[4];
			ContourVertex cv;
			int i,j,ii,jj,ctCrosses,ia,ja;
			for(ii=1;ii<iLen;ii++)
			{
				i = ii-1;
				for(jj=1;jj<jLen;jj++)
				{
					//	In the comments below, let
					//		"Left" refers to lower col indexes of 'Z'
					//		"Right" refers to higher col indexes of 'Z'
					//		"Top" refers to higher row indexes of 'Z'
					//		"Bottom" refers to lower row indexes of 'Z'
					j = jj-1;
					ctCrosses = 0;
					//zmid = 0.25* (Z[ii,jj] + Z[i,jj] + Z[ii,j] + Z[i,j]);
					z00 = Z[ i, j];
					z10 = Z[ii, j];
					z01 = Z[ i,jj];
					z11 = Z[ii,jj];
					//	Left
					crosses[0] = (z00<IsoZ) != (z10<IsoZ);
					//	Right
					crosses[1] = (z01<IsoZ) != (z11<IsoZ);
					//	Top
					crosses[2] = (z10<IsoZ) != (z11<IsoZ);
					//	Bottom
					crosses[3] = (z00<IsoZ) != (z01<IsoZ);
					//	Count crossings
					if(crosses[0]) ctCrosses++;
					if(crosses[1]) ctCrosses++;
					if(crosses[2]) ctCrosses++;
					if(crosses[3]) ctCrosses++;
					//  ctCrosses==0 --> No intersection (do nothing)
					//	ctCrosses==1 --> Not possible (do nothing)
					//	ctCrosses==3 --> Not possible (do nothing)
					//	ctCrosses==4 --> Ambiguous  (do nothing).
					if( ctCrosses==2 )	//	The contour intersects the current cell.
					{
						cv = new ContourVertex(this);
						if( crosses[0] && crosses[1] ) //	Horizontal signature
						{
							//	Left
							dz0 = z10-z00;
							u0 = Math.Abs(z00-IsoZ);
							v0 = Math.Abs(z10-IsoZ);
							u0 = u0/(u0+v0);
							//	Right
							dz1 = z11-z01;
							u1 = Math.Abs(z01-IsoZ);
							v1 = Math.Abs(z11-IsoZ);
							u1 = u1/(u1+v1);
							//	Position
							cv.x = (double)j + 0.5;
							cv.y = (double)i + 0.5*(u0+u1);
							//	Angle
							if(dz0>0.0)
								cv.GradientRadians = Math.Atan2(u1-u0,1.0) + this.RADIANS_1_4;
							else
								cv.GradientRadians = Math.Atan2(u1-u0,1.0) - this.RADIANS_1_4;
							//	Gradient
							cv.dZ_dy = 0.5*(dz0+dz1);
							cv.dZ_dGradient = cv.dZ_dy / Math.Sin(cv.GradientRadians);
							cv.dZ_dx = cv.dZ_dGradient * Math.Cos(cv.GradientRadians);
							//	Storage
							zCrossings[i,j] = cv;
							Vertices.Add(cv);
						}
						else if( crosses[2] && crosses[3] )	//	Vertical signature
						{	
							//	Bottom
							dz0 = z01-z00;
							u0 = Math.Abs(z00-IsoZ);
							v0 = Math.Abs(z01-IsoZ);
							u0 = u0/(u0+v0);
							//	Top
							dz1 = z11-z10;
							u1 = Math.Abs(z10-IsoZ);
							v1 = Math.Abs(z11-IsoZ);
							u1 = u1/(u1+v1);
							//	Position
							cv.x = (double)j + 0.5*(u0+u1);
							cv.y = (double)i + 0.5;
							//	Angle
							if(dz0>0.0)
								cv.GradientRadians = Math.Atan2(u0-u1,1.0);
							else
								cv.GradientRadians = Math.Atan2(u0-u1,1.0) + Math.PI;
							//	Gradient
							cv.dZ_dx = 0.5*(dz0+dz1);
							cv.dZ_dGradient = cv.dZ_dx / Math.Cos(cv.GradientRadians);
							cv.dZ_dy = cv.dZ_dGradient * Math.Sin(cv.GradientRadians);
							//	Storage
							zCrossings[i,j] = cv;
							Vertices.Add(cv);
						}
						else //	Mixed signature
						{
							if( crosses[0] ) // Left
							{
								if( crosses[2] ) // Top
								{
									//	Top-left corner is cut off z10
									//	Top
									dz0 = z11-z10; // dz_dx along top
									u0 = u1 = Math.Abs(z10-IsoZ);
									v0 = Math.Abs(z11-IsoZ);
									u0 = u0/(u0+v0);
									//	Left
									dz1 = z10-z00; // dz_dy along left
									//u1 = Math.Abs(z10-IsoZ);
									v1 = Math.Abs(z00-IsoZ);
									u1 = u1/(u1+v1);
									//	Position  (midpoint of hypotenuse)
									cv.x = (double)j  + 0.5*u0;  // distance from left
									cv.y = (double)ii - 0.5*u1;  // distance from top
									//	Angle
									if( dz1>0.0 )
										cv.GradientRadians = Math.Atan2(u1,u0) + RADIANS_1_4;
									else
										cv.GradientRadians = Math.Atan2(u1,u0) - RADIANS_1_4;
									//	Gradient
									cv.dZ_dx = dz0;
									cv.dZ_dy = dz1;
									if( Math.Abs(dz0) > Math.Abs(dz1) )
										cv.dZ_dGradient = cv.dZ_dx / Math.Cos(cv.GradientRadians);
									else
										cv.dZ_dGradient = cv.dZ_dy / Math.Sin(cv.GradientRadians);
									//	Storage
									zCrossings[i,j] = cv;
									Vertices.Add(cv);
								}
								else if( crosses[3] ) // Bottom
								{
									//	Bottom-left corner is cut off  z00
									//	Bottom
									dz0 = z01-z00; // dz_dx along bottom
									u0 = u1 = Math.Abs(z00-IsoZ);
									v0 = Math.Abs(z01-IsoZ);
									u0 = u0/(u0+v0);
									//	Left
									dz1 = z10-z00; // dz_dy along left
									//u1 = Math.Abs(z00-IsoZ);
									v1 = Math.Abs(z10-IsoZ);
									u1 = u1/(u1+v1);
									//	Position  (midpoint of hypotenuse)
									cv.x = (double)j + 0.5*u0; // distance from left
									cv.y = (double)i + 0.5*u1; // distance from bottom
									//	Angle
									if( dz1>0.0 )
										cv.GradientRadians = RADIANS_1_4-Math.Atan2(u1,u0); //Math.PI-Math.Atan2(u1,u0)-RADIANS_1_4;
									else
										cv.GradientRadians = RADIANS_3_4-Math.Atan2(u1,u0); //Math.PI-Math.Atan2(u1,u0)+RADIANS_1_4;
									//	Gradient
									cv.dZ_dx = dz0;
									cv.dZ_dy = dz1;
									if( Math.Abs(dz0) > Math.Abs(dz1) )
										cv.dZ_dGradient = cv.dZ_dx / Math.Cos(cv.GradientRadians);
									else
										cv.dZ_dGradient = cv.dZ_dy / Math.Sin(cv.GradientRadians);
									//	Storage
									zCrossings[i,j] = cv;
									Vertices.Add(cv);
								}
							}
							if( crosses[1] ) // Right
							{
								if( crosses[2] ) // Top
								{
									//	Top-right corner is cut off z11
									//	Top
									dz0 = z11-z10; // dz_dx along top
									u0 = u1 = Math.Abs(z11-IsoZ);
									v0 = Math.Abs(z10-IsoZ);
									u0 = u0/(u0+v0);
									//	Right
									dz1 = z11-z01; // dz_dy along right
									//u1 = Math.Abs(z11-IsoZ);
									v1 = Math.Abs(z01-IsoZ);
									u1 = u1/(u1+v1);
									//	Position  (midpoint of hypotenuse)
									cv.x = (double)jj - 0.5*u0; // distance from right
									cv.y = (double)ii - 0.5*u1; // distance from top
									//	Angle
									if( dz1>0.0 )
										cv.GradientRadians = RADIANS_1_4-Math.Atan2(u1,u0); //Math.PI-Math.Atan2(u1,u0)-RADIANS_1_4;
									else
										cv.GradientRadians = RADIANS_3_4-Math.Atan2(u1,u0); //Math.PI-Math.Atan2(u1,u0)+RADIANS_1_4;
									//	Gradient
									cv.dZ_dx = dz0;
									cv.dZ_dy = dz1;
									if( Math.Abs(dz0) > Math.Abs(dz1) )
										cv.dZ_dGradient = cv.dZ_dx / Math.Cos(cv.GradientRadians);
									else
										cv.dZ_dGradient = cv.dZ_dy / Math.Sin(cv.GradientRadians);
									//	Storage
									zCrossings[i,j] = cv;
									Vertices.Add(cv);
								}
								else if( crosses[3] ) // Bottom
								{
									//	Bottom-right corner is cut off  z01
									//	Bottom
									dz0 = z01-z00; // dz_dx
									u0 = u1 = Math.Abs(z01-IsoZ);
									v0 = Math.Abs(z00-IsoZ);
									u0 = u0/(u0+v0);
									//	Right
									dz1 = z11-z01; // dz_dy
									//u1 = Math.Abs(z01-IsoZ);
									v1 = Math.Abs(z11-IsoZ);
									u1 = u1/(u1+v1);
									//	Position  (midpoint of hypotenuse)
									cv.x = (double)jj - 0.5*u0;  // distance from right
									cv.y = (double)i  + 0.5*u1;  // distance from bottom
									//	Angle
									if( dz1>0.0 )
										cv.GradientRadians = Math.Atan2(u1,u0) + RADIANS_1_4;
									else
										cv.GradientRadians = Math.Atan2(u1,u0) - RADIANS_1_4;
									//	Gradient
									cv.dZ_dx = dz0;
									cv.dZ_dy = dz1;
									if( Math.Abs(dz0) > Math.Abs(dz1) )
										cv.dZ_dGradient = cv.dZ_dx / Math.Cos(cv.GradientRadians);
									else
										cv.dZ_dGradient = cv.dZ_dy / Math.Sin(cv.GradientRadians);
									//	Storage
									zCrossings[i,j] = cv;
									Vertices.Add(cv);
								}
							}
						}
						//	ContourSegment relations
						ia=i-1;
						ja=j-1;
						if(ia>=0			&& zCrossings[ia, j]!=null)
							Segments.Add(new ContourSegment(this, cv, zCrossings[ia, j]));
						if(ia>=0 && ja>=0	&& zCrossings[ia,ja]!=null)
							Segments.Add(new ContourSegment(this, cv, zCrossings[ia,ja]));
						if(ja>=0			&& zCrossings[ i,ja]!=null)
							Segments.Add(new ContourSegment(this, cv, zCrossings[ i,ja]));
					}
				}
			}
			//	Sort the potential contour segments by ascending affinity
			Segments.Sort();
			//	Segment processing logic.
			i = Segments.Count-1;
			//	Trim any segments that have overlapping vertices
			while(i>=0 && Segments[i].ds<=0.0)
				Segments.RemoveAt(i--);
			//	Join the unjoined vertices.
			while(i>=0)
				Segments[i--].JoinUnjoinedVertices();
			//	Contour formation logic, begin by sorting vertices by descending LowestAffinity.
			Vertices.Sort();
			Vertices.Reverse();
			int nMax = Vertices.Count;
			j=0;
			for(i=0; i<Vertices.Count && nMax>3; i++)
			{
				if(Vertices[i].IdContour==int.MinValue)
				{
					Contours.Add(new Contour(j, Vertices[i], nMax, MeshRect));
					nMax -= Contours[j].Vertices.Length;
					j++;
				}
			}
			//	Sort by contour length descending (longest contour first, shortest last)
			Contours.Sort();
			Contours.Reverse();
			//	Adjust the vertex coordinates to the DataRect
			for(i=0; i<Contours.Count; i++)
				Contours[i].IdContour = i;
		}
		private double comparisonZ;
		int IComparable.CompareTo(object obj)
		{
			comparisonZ = ((IsoContours)obj).IsoZ;
			if(IsoZ > comparisonZ)
				return 1;
			if(IsoZ == comparisonZ)
				return 0;
			return -1;
		}
		/// <summary>
		/// The pen used for drawing the axis, width is automatically scaled by DrawScale.
		/// </summary>
		public Pen Pen
		{
			get { return pen; }
			set
			{
				pen = value;
				scaledPen = new Pen(pen.Color,pen.Width*drawScale);
			}
		}
		private Pen pen = new Pen(Color.Black,1.0f);
		private Pen scaledPen = new Pen(Color.Black,1.0f);
		/// <summary>
		/// The scale with which to draw the axis.  Automatically scales the font and the pen.
		/// </summary>
		public float DrawScale
		{
			get { return drawScale; }
			set
			{
				if (value!=drawScale)
				{
					drawScale = value;
					scaledPen = new Pen(pen.Color,pen.Width*drawScale);
				}
			}
		}
		private float drawScale = 1.0f;
	}
	/// <summary>
	/// Represents the affinity with which two points (A and B) can join to form a contour segment.
	/// </summary>
	public class ContourSegment : IComparable
	{
		/// <summary>
		/// The parent of this object.
		/// </summary>
		public IsoContours Parent;
		/// <summary>
		/// The starting vertex.
		/// </summary>
		public ContourVertex A;
		/// <summary>
		/// The ending vertex.
		/// </summary>
		public ContourVertex B;
		/// <summary>
		/// A positive number in [0, inf).  The higher this number, the more likely that two adjacent
		/// vertices should be grouped.
		/// The affinity is a function of the agreement in angles and the distance between vertices.
		/// It's just a formula that I made up ad-hoc.  It should perform pretty well, but there may
		/// be room for improvement.
		/// </summary>
		public double Affinity = Double.NegativeInfinity;
		/// <summary>
		/// The change in x moving from A to B.
		/// </summary>
		public double dx;
		/// <summary>
		/// The change in y moving from A to B.
		/// </summary>
		public double dy;
		/// <summary>
		/// The distance between vertices A and B.
		/// </summary>
		public double ds;
		/// <summary>
		/// The change in x per unit distance moving from A to B.
		/// </summary>
		public double dx_ds;
		/// <summary>
		/// The change in y per unit distance moving from A to B.
		/// </summary>
		public double dy_ds;
		/// <summary>
		/// The change in Z per unit of the vector A-->B rotated by +90 degrees;  This is calculated
		/// by rotating this.dZ_dx and this.dZ_dy relative to A-->B such that A-->B is at 0 degrees.
		/// </summary>
		public double dZ_dCosAB, dZ_dSinAB;
		/// <summary>
		/// This is calculated by averaging dZ_dx and dZ_dy from A and B.
		/// </summary>
		public double dZ_dx, dZ_dy;
		/// <summary>
		/// The contour segment joining vertices A and B moving from A to B.  The constructed segment will
		/// map the arguemnt A to either this.A or this.B and map the argument B to the other member
		/// this.A or this.B such that dZ_dSinAB less or equal 0.0.  Thus, the segment A-->B is always 
		/// compatible with a clockwise polygon (where Z is larger inside the polygon) and where A.Next is B
		/// and B.Prev is A.
		/// </summary>
		/// <param name="A">The starting vertex.</param>
		/// <param name="B">The ending vertex.</param>
		public ContourSegment(IsoContours Parent, ContourVertex A, ContourVertex B)
		{
			this.Parent = Parent;
			dx = B.x-A.x;
			dy = B.y-A.y;
			ds = Math.Sqrt(dx*dx+dy*dy);
			dx_ds = dx/ds;
			dy_ds = dy/ds;
			dZ_dx = 0.5*(A.dZ_dx+B.dZ_dx);
			dZ_dy = 0.5*(A.dZ_dy+B.dZ_dy);
			dZ_dCosAB =  dZ_dx*dx_ds + dZ_dy*dy_ds;
			dZ_dSinAB = -dZ_dx*dy_ds + dZ_dy*dx_ds;
			if(dZ_dSinAB<=0.0)
			{
				this.A = A;
				this.B = B;
			}
			else
			{
				this.A = B;
				this.B = A;
				dx = -dx;
				dy = -dy;
				dx_ds = -dx_ds;
				dy_ds = -dy_ds;
				dZ_dx = -dZ_dx;
				dZ_dy = -dZ_dy;
				dZ_dCosAB = -dZ_dCosAB;
				dZ_dSinAB = -dZ_dSinAB;
			}
			this.Affinity = (1.0+(A.dZ_dx*B.dZ_dx+A.dZ_dy*B.dZ_dy)/(A.dZ_dGradient*B.dZ_dGradient))/ds;
			//this.Affinity = 1.0/ds;
		}
		/// <summary>
		/// Join vertices by linking them in a double-linked-list and recording the linking affinity.
		/// </summary>
		public void JoinVertices()
		{
			this.A.Next = this.B;
			if(this.Affinity > (double)this.A.HighestAffinity)
			{
				this.A.LowestAffinity = this.A.HighestAffinity;
				this.A.HighestAffinity = (float)this.Affinity;
			}
			else
			{
				if(this.Affinity > (double)this.A.LowestAffinity)
					this.A.LowestAffinity = (float)this.Affinity;
			}

			this.B.Prev = this.A;
			if(this.Affinity > (double)this.B.HighestAffinity)
			{
				this.B.LowestAffinity = this.B.HighestAffinity;
				this.B.HighestAffinity = (float)this.Affinity;
			}
			else
			{
				if(this.Affinity > (double)this.B.LowestAffinity)
					this.B.LowestAffinity = (float)this.Affinity;
			}
		}
		/// <summary>
		/// Join vertices if they are currently unjoined.
		/// </summary>
		public void JoinUnjoinedVertices()
		{
			if(this.A.Next==null && this.B.Prev==null)
				this.JoinVertices();
		}
		private double comparisonAffinity;
		int IComparable.CompareTo(object obj)
		{
			comparisonAffinity = ((ContourSegment)obj).Affinity;
			if(Affinity > comparisonAffinity)
				return 1;
			if(Affinity == comparisonAffinity)
				return 0;
			return -1;
		}
		public PointF[] Paintable(PaintingInfo pi)
		{
			PointF[] output = new PointF[2];
			float dwi = pi.Rect.Width/Parent.MeshRect.Width;
			float dhi = pi.Rect.Height/Parent.MeshRect.Height;
			output[0].X =				 ((float)A.x-Parent.MeshRect.X)*dwi + pi.Rect.X;
			output[0].Y = pi.Rect.Height-((float)A.y-Parent.MeshRect.Y)*dhi + pi.Rect.Y; // Y is inverted
			output[1].X =				 ((float)A.x-Parent.MeshRect.X)*dwi + pi.Rect.X;
			output[1].Y = pi.Rect.Height-((float)A.y-Parent.MeshRect.Y)*dhi + pi.Rect.Y; // Y is inverted
			return output;
		}
	}
	/// <summary>
	/// The vertex of a contour.
	/// </summary>
	public class ContourVertex : IComparable
	{
		/// <summary>
		/// The parent of this object.
		/// </summary>
		public IsoContours Parent;
		/// <summary>
		/// Horizontal zero-based floating point number that floats between the zero-based integer indexes of Z.
		/// </summary>
		public double x;
		/// <summary>
		/// Vertical zero-based floating point number that floats between the zero-based integer indexes of Z.
		/// </summary>
		public double y;
		/// <summary>
		/// The direction towards which z increases (in radians).
		/// </summary>
		public double GradientRadians;
		/// <summary>
		/// The gradient magnitude (this is positive only because the gradient direction is pointed towards the positive gradient).
		/// </summary>
		public double dZ_dGradient;
		/// <summary>
		/// The local change in Z as x increases.
		/// </summary>
		public double dZ_dx;
		/// <summary>
		/// The local change in Z as y increases.
		/// </summary>
		public double dZ_dy;
		/// <summary>
		/// The ID number of the contour to which this vertex belongs.
		/// </summary>
		public int IdContour = int.MinValue;
		/// <summary>
		/// The previous contour vertex.  As the linked list is traversed from Prev to Next,
		/// Z increases in the clockwise direciton and decreases in the counter-clockwise direction.
		/// </summary>
		public ContourVertex Prev;
		/// <summary>
		/// The next contour vertex.
		/// </summary>
		public ContourVertex Next;
		/// <summary>
		/// The higher of two affinities.
		/// </summary>
		public float HighestAffinity = float.MinValue;
		/// <summary>
		/// The lower of two affinities.
		/// </summary>
		public float LowestAffinity = float.MinValue;
		public ContourVertex(IsoContours Parent)
		{
			this.Parent = Parent;
		}
		/// <summary>
		/// A string representation of this thing.
		/// </summary>
		/// <returns></returns>
		public string ToString()
		{
			return "["+
				"("+ ((float)x).ToString() +","+ ((float)y).ToString() + ")"
				+" , "+
				"("+ ((float)dZ_dx).ToString() +","+ ((float)dZ_dy).ToString() + ")"
				+"]";
		}
		private float CompareTo_LowestAffinity;
		/// <summary>
		/// Sort by the lowest affinity.
		/// </summary>
		/// <param name="o"></param>
		/// <returns></returns>
		int IComparable.CompareTo(object o)
		{
			CompareTo_LowestAffinity = ((ContourVertex)o).LowestAffinity;
			if(this.LowestAffinity>CompareTo_LowestAffinity)
				return 1;
			if(this.LowestAffinity==CompareTo_LowestAffinity)
				return 0;
			return -1;
		}
		public PointF[] Paintable(PaintingInfo pi)
		{
			PointF[] output = new PointF[2];
			float dwi = pi.Rect.Width/Parent.MeshRect.Width;
			float dhi = pi.Rect.Height/Parent.MeshRect.Height;
			output[0].X =				 ((float)this.x-Parent.MeshRect.X)*dwi + pi.Rect.X;
			output[0].Y = pi.Rect.Height-((float)this.y-Parent.MeshRect.Y)*dhi + pi.Rect.Y; // Y is inverted
			output[1].X =				 ((float)this.x-Parent.MeshRect.X+(float)Math.Cos(GradientRadians)*5.0f)*dwi + pi.Rect.X;
			output[1].Y = pi.Rect.Height-((float)this.y-Parent.MeshRect.Y+(float)Math.Sin(GradientRadians)*5.0f)*dhi + pi.Rect.Y; // Y is inverted
			return output;
		}
	}
	/// <summary>
	/// A contour.
	/// </summary>
	public class Contour : IComparable
	{
		/// <summary>
		/// The ID tag associated with this contour.  When you set this property, the IdContour property
		/// is automatically set for all elements in this.Vertices.
		/// </summary>
		public int IdContour
		{
			get
			{
				return this.idContour;
			}
			set
			{
				this.idContour = value;
				for(int i=0; i<Vertices.Length; i++)
					Vertices[i].IdContour = value;
			}
		}
		private int idContour = int.MinValue;
		public bool IsClosed = false;
		public RectangleF DataRect;
		public ContourVertex[] Vertices;
		/// <summary>
		/// Construct a new contour with the specified ID 
		/// </summary>
		/// <param name="IdContour">An integer ID number associated with this contour.</param>
		/// <param name="cv">A ContourVertex in the contour.</param>
		/// <param name="nMax">The maximum number of vertices that the contour could contain.</param>
		public Contour(int IdContour, ContourVertex cv, int nMax, RectangleF DataRect)
		{
			this.DataRect = DataRect;
			this.idContour = IdContour;
			if(cv.IdContour>=0 && cv.IdContour!=idContour)
				Vertices = new ContourVertex[0];
			int ct = 1;
			ContourVertex cvFirst = cv;
			ContourVertex cvLast = cv;
			//	Make cvLast keep going until end.
			while(cvLast.Next!=null)
			{
				if(cvLast.Next == cvFirst)
				{
					this.IsClosed = true;
					break;
				}
				else
				{
					ct++;
					cvLast = cvLast.Next;
				}
				if(ct>nMax)
					throw new Exception("Loop continues too far.");
			}
			if(!IsClosed)
			{
				while(cvFirst.Prev!=null)
				{
					ct++;
					cvFirst = cvFirst.Prev;
					if(ct>nMax)
						throw new Exception("Loop continues too far.");
				}
			}
			Vertices = new ContourVertex[ct];
			cv = cvFirst;
			for(int i=0; i<ct; i++)
			{
				Vertices[i] = cv;
				cv.IdContour = IdContour;
				cv = cv.Next;
			}
		}
		private PointF[] paintCache;
		private RectangleF lastPaintRect = RectangleF.Empty;
		public PointF[] Paintable(PaintingInfo pi)
		{
			if(paintCache==null || paintCache.Length!=Vertices.Length)
			{
				lastPaintRect = RectangleF.Empty;
				paintCache = new PointF[Vertices.Length];
			}
			if( !lastPaintRect.Equals(pi.Rect) )
			{
				lastPaintRect = pi.Rect;
				PointF pf = PointF.Empty;
				ContourVertex cv;
				float dwi = pi.Rect.Width/DataRect.Width;
				float dhi = pi.Rect.Height/DataRect.Height;
				for(int i=0; i<Vertices.Length; i++)
				{
					cv = Vertices[i];
					pf.X =				  ((float)cv.x-DataRect.X)*dwi + pi.Rect.X;
					pf.Y = pi.Rect.Height-((float)cv.y-DataRect.Y)*dhi + pi.Rect.Y; // Y is inverted
					paintCache[i] = pf;
				}
			}
			return paintCache;
		}
		private int comparisonLength;
		int IComparable.CompareTo(object obj)
		{
			comparisonLength = ((Contour)obj).Vertices.Length;
			if(Vertices.Length > comparisonLength)
				return 1;
			if(Vertices.Length == comparisonLength)
				return 0;
			return -1;
		}
	}
}